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Abstract. - We study analytically the behavior of the largest Lyapunov exponent Ai for a one- 
dimensional chain of coupled nonlinear oscillators, by combining the transfer integral method 
and a Riemannian geometry approach. We apply the results to a simple model, proposed for 
the DNA denaturation, which emphasizes a first order-like or second order phase transition 
depending on the ratio of two length scales: this is an excellent model to characterize Ai as a 
dynamical indicator close to a phase transition. 



Introduction. - Phase transition plays a central role in equilibrium and non equilibrium 
statistical physics books and lectures in particular because it exemplifies the paradigmatic 
concept of universality in physics. In the theory of dynamical systems, the concept of Lyapunov 
exponent has also attracted a lot of attention because it defines unambiguously a sufficient 
condition for chaotic instability, but unfortunately except for very few systems it is already 
an extremely difficult task to derive analytically the expression of the largest one, Ax, as a 
function of the energy density. As some promising results have been recently obtained to 
describe some properties of high-dimensional dynamical systems j3|J(|, by combining tools 
developed in the framework of dynamical systems with concepts and methods of equilibrium 
I 1 statistical mechanics, the idea that both concepts could be related was proposed recently 0. 

During the last years, the process by which two strands of DNA unbind upon heating, 
called DNA denaturation or melting, has motivated a lot of works |]8|-|ll| and in particular an 
\ extremely simplified dynamical model ||l2f was proposed and studied. This one-dimensional 
L" ' hamiltonian has the following expression 

X ; H = £ ^yl + | (y n - y n+l f + D ( e -* - l) 2 (1) 

where m corresponds to the effective mass of nucleotides, K the coupling constant and D (resp. 
a) the depth (resp. inverse length scale) of the Morse potential which mimic the interactions 
between groups of atoms of opposite strands. 
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This potential with only nearest-neighbor interactions is well suited for using the transfer 
integral method in order to derive the canonical partition function. It can be shown JT^ | 
that this system can be mapped to the quantum mechanical analogy of a particle in a Morse 
potential and that all thermodynamical quantities could finally be expressed as functions of 
the eigenvalues e q and eigenfunctions 4> q of this Schrodinger problem. In particular, we showed 
the existence of a critical temperature T c — 2\/2KD/aks corresponding to a second order 
phase transition between the so-called native (or double strand) state with the particles in the 
bottom of the Morse well, and the denaturated state, with the particles on the Morse plateau. 

The dynamical properties of this one-dimensional model are also very interesting and have, 
in particular, emphasized the role of localized oscillating excitations, called discrete breathers, 
as precursors effects driving this phase transition. Spatiotemporal studies of the dynamics 
reveals also intermittency-like features and has led us to consider the importance of its chaotic 
properties as an important ingredient to characterize and to explain this dynamical instability. 
This is the reason why it is important to study the Lyapunov behavior as a function of the 
energy density, not only by computing it numerically but also (if possible !) by deriving its 
expression using the Riemannian geometry approach proposed recently . 

Riemannian Geometry Approach. - The main idea is that the chaotic hypothesis is at 
the origin of the validity of equilibrium statistical physics, and this should be traced somehow 
in the dynamics and therefore in the largest Lyapunov exponent. Through a reformulation 
of Hamiltonian dynamics in the language of Riemannian geometry Q , the method proposes 
therefore to relate the microscopic dynamics to the statistical averages fL3f . Applied to a 
typical high-dimensional hamiltonian system, this corresponds to study the parametric-like 
instability of geodesies, driven by the fluctuations of the curvature. Indeed, chaos can be 
induced not only by negative curvatures but also by positive ones, provided they are fluctu- 
ating, as it is less well known. Once the curvature kq and its fluctuations a K are accurately 
determined, the method uses a gaussian stochastic process to approximate the effective curva- 
ture felt along a geodesic and finally one ends up with the following expression of the largest 
Lyapunov exponent 




2 V 3A 

In this definition, r, the relevant time scale associated to the stochastic process is function 
of the two following timescales: t\ ~ 7t/2^/ko + cr K is the time needed to cover the distance 
between two successive conjugate points along the geodesies, whereas T2 ~ y/Ko/cr K is related 
to the local curvature fluctuations. The general rough physical estimate t~(1/ti + 1/t2) 1 
completes finally the analytical estimate of Ai that we would like to continue now by the 
calculation of the mean value of the curvature and its fluctuations as a function of the energy 
density. 

Curvature and fluctuations along geodesies. Using the simplifying Eiscnhart metric, 
the Ricci curvature Kjj(y) corresponds to the Laplacian of the potential energy. Defining the 
function g(y) — 2e~ 2ay — e~ ay , we get the following expression 

K ° = A^f" ~ 2K + 2a2D ^(v))^ = 2K+ 2a 2 D (l - if T < T c , (3) 

= 2K if T>T C (4) 
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Fig. 1 - Evolution of ko (circles) and a K (squares) as a function of the temperature T. Circles and 
squares correspond to microcanonical results, the solid lines to results obtained with the transfer 
integral method; the dashed line to the analytical expression (^). Finally, the dotted line correspond 
to the low temperature approximation (J12[) of the fluctuations of curvature. 



where we have used the equality of mean values in microcanonical and canonical («) ca „ 
ensembles | fl5f . The last expression is in particular derived using the expression of the ground- 
state 4>o{y) of the quantum mechanical analogy since (g(y)) C an = (</ , o|5(j/)|0o)can- Figure [l] 
attests that the agreement is excellent between the above analytical expression, the canonical 
transfer integral results and microcanonical molecular dynamics simulations, obtained using 
the best 4th order symplectic integrator due to MacLahan-Attela j^J. 

It is interesting to note that, the mean curvature is always positive whereas the local one 
is negative close to the inflexion point of the Morse potential. In addition, expression (||) 
emphasizes the role of discreteness since the curvature would be almost constant in the con- 
tinuum limit a 2 D <C K. This will have further consequences on the evolution of Ai when the 
discreteness is changed. 

Fluctuations of curvature. Since numerical computations are simpler in the micro- 
canonical ensemble, while analytical calculations are simpler in the canonical one, to get the 
fluctuations of the curvature in the microcanonical ensemble, one needs to add the following 
corrective term ]l5| 

= (S 2 K R ) can - ( ^#^V . (5) 

The expression of the energy 



dp J V 9(3 



U 1 1 dhiZ cm T 2 

— = t- 2 ^- = k B T + D — if T <T C , 6 

N 2(3 N d(3 T 2 ' w 

= k B T if T>T C , (7) 
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emphasizes the second order of the phase transition, and finally using the expression of the 
mean value of the curvature, we obtain the microcanonical/canonical correction 

4« 4 D 2 /V T 2 

(S 2 K R ), (6 2 K R ) can = - " t= ( 8 ) 

1T V 2K T c 

valid below T c . This quadratic correction vanishes for very low temperature and converges 
toward a negative constant close to T c . 

The expression of the fluctuations of curvature in the canonical ensemble is then 

(S 2 K R ) can = (K 2 H (y) - (K R (y)) 2 ) can = Aa 4 D 2 ^ ({g{y l )g{y 3 )) can ~ (<?(</;)>can(<?(%)>c« 

The second part of the parenthesis could be easily computed and gives a term N 2 (g} 2 . The 
first one can be evaluated with the transfer integral method since, thanks to the periodic 
boundary condition, one needs to take into account only the difference in indices (i — j). 
However, one cannot integrate directly the integral over the resulting (N — 1) variables, and 
one needs to decompose twice on the orthonormal basis {4> q } of the transfer operator. Defining 
the matrix elements of g(y) in this base 

/+oo 
r qo (y)g(y)<t>k{y)dy = (K\9(y)\<t>k) , (9) 
-oo 

we finally obtain 

(9(y P )9(y N )) ca n = L^e-^-^-^ e~ N ^ |M M | 2 . 
c i 

K,q 



As only the lowest eigenvalue Sq will contribute in the thermodynamic limit, we finally get 

_|Mo, ? | 2 
1 - 

9=1 



(5 2 K R)can = 4a^ E J^l_ £o) (10) 



by recognizing a geometrical sum. Figure (|l|) emphasizes the excellent agreement between 
microcanonical simulations (squares) and the above analytical expression (solid line) where 
we have used the eigenfunctions of the transfer integral operator. 

Low temperature approximation. The above expression ( |Io|) is however difficult to 
handle in general but could be simplified in the low temperature regime since the gap between 
the eigenvalue of the ground state eo and the other ones justify to neglect the exponential term 
in the denominator. This approximation breaks down of course close to the phase transition 
since the critical temperature is defined by the disappearance of the ground state eo in the 
continuum. We obtain thus 

{5 2 K R ) can ^ Aa i D 2 NY J \{^W q )\ 2 =^D 2 N[{g 2 )-{g) 2 ] . (11) 

■7=1 

Using the procedure described for the calculation of (g), one obtains the formula (g 2 ) = 
(l — ^1 + ^ + ^5-)j where the prefactors 2 and 4 have been renormalized to be in 
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agreement with numerical simulations; the discreteness of the chain, which drives the phase 
transition, is the main reason for this renormalization. Adding finally the microcanonical 
correction (fi), we obtain 



(S 2 K f 



N - 1 



4a 4 L> 



T 



5T 



2T 2 
If 




(12) 



The dotted line plotted in Fig. (|l|) attests the quality of this low temperature approximation 
which leads to a linear dependence versus temperature of the fluctuations a 2 ~ 20a 4 D 2 T/T c 
at the lowest order. 

In the neighborhood of the phase transition, the last bound state is too close to the 
continuum: the exponential term in the denominator of Eq. (|To|) cannot be neglected and 
one needs to take into account fully the interactions with the continuum. Due to the rapid 
oscillations of the scattering eigenstates, one could use the approximate expression for the 
phase shift of the scattering states, as it was successfully shown in the calculation of the static 
structure factor and therefore for the susceptibility as its limit Jl^| . But, unfortunately for 
the present calculation, there is an important contribution of eigenstates with nonvanishing 
eigenvectors which prevents this high temperature approximation to succeed. 

Largest Lyapunov exponent. Once the curvature (||) and its fluctuations jic| ) have 
been calculated, it is straightforward to compute the two timescale T\ and T2- Having in 
mind that the estimate of the decorrelation time scale r is still somehow rudimentary j?J in 
the Riemannian framework outlined above, Fig. (|^) shows that the alternative definition r ~ 
(1/ti + 5/T2) 1 gives a particularly striking agreement between the microcanonical numerical 
results obtained using the standard algorithm Jl6| and the analytical derivation proposed 
here. Let us emphasize that the agreement is excellent on the whole interval contrary to 
earlier results obtained for the 4 chain (jl7|. The approach described here has a domain of 
validity even larger than expected since the agreement is very good even in the region where 
the fluctuations of curvature are of the order or greater than the curvature itself. 

It is also very interesting to remark that the simultaneous use of the Riemannian geometry 
approach and the transfer integral method, that we propose here, could be easily extended to 
the treatment of any chain of oscillators with on-site potential and nearest-neighbor harmonic 
coupling. Once the transfer integral operator is solved analytically as it is possible here for 
Hamiltonian (|l|) or numerically in general, the evolution of the largest maximal Lyapunov 
exponent as a function of the energy density is directly derived from expressions of the curva- 
ture {^) and fluctuations (10), by simply replacing 2a 2 Dg(y) by the second derivative of the 
onsite potential. 

A closer look at very low energy density, shows however a non surprising disagreement 
since the inset plotted in Fig. (^J) clearly violates the low energy density approximate law 
Ai ~ &k> which would give here a linearly increasing function of the temperature according 
to the lowest order estimate of a K . The reason is presumably the breakdown of the stochastic 
approximation for the average curvature in this very low energy limit where the system is 
almost quadratic and therefore the crucial ergodic hypothesis is not completely fullfilled any 
more. More interestingly, the numerical model emphasizes an unusual 3/2 power law contrary 
to the generally reported u 2 dependence of Ai in high-dimensional dynamical systems |],|t]]. 

Let us note that the expression of kq is strongly dependent of the discreteness parameter 
a 2 D/K, but the fluctuations a K are not. As Fig. (0) shows that, away from the critical region, 
smaller the curvature is, greater the Lyapunov is, one can expect that the more discrete the 
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Fig. 2 - Evolution of the maximal Lyapunov exponent as a function of the energy density u = U /TV. 
The symbols correspond to microcanonical numerical simulations for TV =500 (triangles), TV = 10 3 
(squares), TV = 2.10 3 (crosses) and TV = 10 4 (diamonds). The solid line corresponds to the analytical 
estimation using the transfer integral method whereas the dashed line corresponds to the analytical 
expression valid at low energy. The vertical dotted line shows the position of the critical energy 
density u c — ksT c + D. In the inset, we plot Ai as a function of the energy density u in log-log scale. 
The symbols correspond to microcanonical numerical simulations for TV = 500 whereas the solid line 
to a 3/2 power law. 

Fig. 3 - Lyapunov exponent and order of the phase transition. The solid line shows the evolution 
of the harmonic model with second order phase transition, whereas the symbols are referring to the 
model with anharmonic coupling with a first-order like phase transition. The number of sites in the 
chain is TV =500 (crosses), TV = 10 3 (diamonds) and TV = 10 4 (circles). The vertical dotted (resp. 
dash-dotted) line corresponds to the critical energy density of Hamiltonian with harmonic (resp. 
anharmonic) coupling. 



chain is, the more chaotic its dynamics would be. This result is clearly in qualitative agree- 
ment with the discovery that spontaneously created discrete breathers are actually chaotic 
excitations and that their domain of stability is far from the continuum limit fl§| ]. 

Lyapunov and Phase transition. - This model is also a very interesting case where one 
can characterize Ai as a dynamical indicator of a phase transition. In this case of a second 
order phase transition, one obtains a critical slowing down reminiscent of the results obtained 
for the Heisenberg Mean-field model derived by Firpo ||. We would like also to emphasize 
the singular behavior of kq, a K and Ai at the critical energy density, in complete agreement 
with the exciting conjectures linking them to a topology change in the underlying manifold, 
being itself an indicator of a thermodynamic phase transition. 

More interestingly, we also have studied Hamiltonian ([!]) whith an anharmonic coupling 
potential y[l + e~ a ^ Vn+yn - 1 ' > ] (y n — y n -i) 2 which has the property to describe the varying 
backbone stiffness of the DNA. This hamiltonian was shown to have a first order like phase 
transition jl9| when the ratio of the two inverse length scales a/ a is lower than 1/2. If 
analytical estimate of the Lyapunov exponent are not possible, Fig. |^ shows its evolution 
obtained using microcanonical molecular dynamics simulations. In the low energy limit, the 
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additional parenthesis do not produce of course modifications, but it emphasizes an abrupt 
change pp| of the Lyapunov exponent at the critical energy density as if Ai was a dynamical 
order parameter indicating the first order phase transition. 

Conclusion. - In this letter, we have presented one of the very few analytical calculation 
of the largest Lyapunov exponent in a high-dimensional dynamical system p|. As expected, 
in addition to serve as a criterion for chaos, or as a characteristic time scale of chaoticity, Ai is 
an excellent dynamical indicator of the presence of the phase transition and could be used not 
only to describe the dynamics but also the statistical properties of high dimensional systems. 

It is also interesting to consider the behavior of another important dynamical indicator, 
the participation ratio of the normalized Lyapunov vector V\ associated to the maximal Lya- 
punov Ai. Defined as £ = 1/ (NJ2iLi [^i (*) 2 +V\(i + N) 2 ] , where the first (resp. last) N 
components are associated to the evolution of linear perturbations of y n (resp. y n ) in tangent 
space, £ is an indicator of localization (2l], : it is of order one if the vector is extended and 
of order 1/N if localized. Here, we found that the phase transition corresponds to a crossover 
from a localized state in tangent space at low energy density to a more extended state just after 
the phase transition, confirming that the tangent space trough the largest Lyapunov exponent 
Ai, but also its associated eigenvector Vi through the participation ratio £, are surprisingly 
good dynamical indicators for emphasizing thermodynamical phase transitions. 

* * * 

We wish to acknowledge very helpful discussions with A. Alastuey, J. Farago, M.-C. Firpo, 
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